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Close-packed, classical dimer models on three-dimensional, bipartite lattices harbor a Coulomb phase with 
power-law correlations at infinite temperature. Here, we discuss the nature of the thermal phase transition 
out of this Coulomb phase for a variety of dimer models which energetically favor crystalline dimer states 
with columnar ordering. For a family of these models we find a direct thermal transition from the Coulomb 
phase to the dimer crystal. While some systems exhibit (strong) first-order transitions in correspondence with 
the Landau-Ginzburg- Wilson paradigm, we also find clear numerical evidence for continuous transitions. A 
second family of models undergoes two consecutive thermal transitions with an intermediate paramagnetic 
phase separating the Coulomb phase from the dimer crystal. We can describe all of these phase transitions in 
one unifying framework of candidate field theories with two complex Ginzburg-Landau fields coupled to a U(l) 
gauge field. We derive the symmetry-mandated Ginzburg-Landau actions in these field variables for the various 
dimer models and discuss implications for their respective phase transitions. 



PACS numbers: 05.30.-d, 02.70.Ss, 64.60.-i, 71.10.Hf 
I. INTRODUCTION 

Constraints are a pervasive feature of strongly correlated 
systems. For instance, in Mott insulators, a large atomic 
Coulomb repulsion effectively constrains the charge of each 
ion to be fixed, while still allowing spin and orbital fluctu- 
ations. This situation, in which the dominant terms in the 
Hamiltonian impose constraints, but local fluctuations remain 
strong, provides a challenge to physical understanding. In 
frustrated magnets, it is common to observe a "cooperative 
paramagnetic" regime, in which the dominant exchange inter- 
actions impose strong constraints on the spin configurations, 
but the spins still manage to remain strongly fluctuating. The 
"spin ice" materials, in which rare earth Ising moments locally 
satisfy Pauling's ice rules, provide a prominent and beautiful 
set of examples, which have stimulated a rich interplay be- 
tween theory and experiment. Constraints on the spin phase 
space have been implicated in the physics of diverse other 
magnetic materials, such as the spinel chromite and the A- 
site diamond antiferromagnetic spinelJ^. 

Generally, residual interactions, subdominant to those re- 
sponsible for the constraints, lead to a quenching of the re- 
maining fluctuations. To quantify this, we may associate the 
dominant interactions with a temperature, Tq, below which 
the constraints are well satisfied and the system is highly cor- 



related. We will assume that fluctuations amongst the con- 
strained states are removed at another temperature Tc <C To, 
which is determined by subdominant effects. Often this 
quenching of the constrained fluctuations is associated with 
a symmetry breaking, such as magnetic ordering or lattice 
deformation. This phase transition occurs in a very dif- 
ferent environment from conventional order-disorder tran- 
sitions, in which the high temperature phase is a weakly 
coiTelated paramagnet. Here, the strong constraints imply 
strong correlations in the cooperative paramagnet. It has re- 
cently been appreciated th at suc h coiTelations can drastically 
affect phase transition(s).^"^^^ Transitions can be induced 
where none would otherwise be present, and furthermore, 
symmetry-mandated transitions may be modified from their 
usual Landau-Ginzburg-Wilson (LGW) universality classes. 

In this paper, we explore these phenomena in a large set 
of classical dimer models on the cubic lattice. Such dimer 
models are defined by a constrained phase space consisting of 
close-packed dimer coverings, in which dimers occupy (some) 
links of the lattice, and the constraint is that each site is over- 
lapped by one and only one dimer. In these models, the con- 
straint is exactly satisfied, corresponding to the limit Tq oo. 
We expect that more realistic models are well approximated 
by this situation provided Tc ^Tq. For models with a gap (of 
order 0{kBTi^)) to states violating the constraint (as in spin 
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ice), the approximation is in fact exponentially good, since vi- 
olations of the constraint occur with an Arrhenius probability 
a exp(-ro/Tc). 

In the cubic dimer models we study, a great deal is under- 
stood about the nature of the constraint induced correlation^^] 
which have a power-law form. This can be cast (see Sec. |IV| i 
into a sort of pseudo-dipolar form, leading to the name 
"Coulomb phase" for the high temperature T > Tc region. 
Moreover, the dipolar correlations can be identified with an 
emergent Coulomb gauge field, similar to that appearing in 
true electromagnetism. Coulomb phases arise in a variety of 
other context^^, for instance, the ice rules constraint related 
to spin ice also leads to a Coulomb phase^-i^" For these 
types of systems, the gauge description has lead to some the- 
oretical progress in understanding the consequent unconven- 
tional criticality. In such cases, the transitions are expected to 
involve dual "monopole" fields, which couple to the emergent 
gauge field and carry the associated gauge charge. The full 
field theory therefore has a multicomponent Ginzburg-Landau 
form. The monopoles are "fractional" degrees of freedom 
in the sense that the symmetry-breaking order parameters (if 
any) are composites of these fields. We outline the derivation 
of this result in Sec.lIVI 

Though this construction of non-LGW critical theories has 
been discussed in some isolated instances previously, an un- 
equivocal verification of the theory has proved difficult. In 
particular, the numerical experiments in Refs. ,4J3. while pro- 
viding evidence for unconventional criticality, are not in good 
quantitative agreement with the theoretical expectations. In 
particular, the numerical estimates'* of various critical expo- 
nents, e.g. I' = 0.50(4), a = 0.56(7) and 77 = -0.02(5), are 
indicative of an unexplained tricritical behavior. 

In this paper, we perform a much more systematic investi- 
gation of a range of dimer models, in order to test the theo- 
retical picture on a grander scale, in which many qualitative 
comparisons are possible. We find that the gauge theory does 
an excellent job on these qualitative tests, providing under- 
standing of the numerical results for nearly all cases. 



A. Outline of models and results 

Before going into a detailed discussion we first provide a 
brief overview of the models we will study and our main re- 
sults. We start by introducing a family of close-packed, classi- 
cal dimer models on the cubic lattice. The elementary degrees 
of freedom in these models are hard-core dimers occupying 
the bonds of the cubic lattice with the constraint that every 
site in the lattice is part of exactly one such dimer We further 
introduce a (potential) energy scale favoring dimer coverings 
with parallel dimers on neighboring bonds. At low temper- 
ature these models exhibit long range columnar order of the 
dimers. In general, there are six distinct columnar ordering 
patterns that can maximize the number of parallel dimers on 
neighboring bonds as illustrated in Fig. [T] 

The main distinction between the various models we will 
consider then comes from a selection of a subset of these en- 
ergetically favored ordering patterns. The Hamiltonian of a 
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FIG. 1 : (color online) The six dimer coverings with columnar dimer 
ordering that maximize the number of parallel dimers on neighbor- 
ing bonds of the cubic lattice. Our family of classical dimer models 
energetically favors these ordering patterns with the '6-GS' model 
favoring all six states, the '4-GS' model favoring states 3 — 6, the 
'2-GS' model favoring states 1 — 2 and the '1-GS' model favoring a 
single state only. Further variations are discussed in the text. 



parent model that favors all six possible ordering patterns and 
which we therefore call the '6-GS' model is given by 



n 



6-GS 



V/ 



+ n||). 



(1) 



□ 



Here n//, and n\ \ count the number of parallel dimers on 
neighboring bonds along the x^y,z lattice directions and the 
sum runs over all square plaquettes of the cubic lattice. 

Our first family of dimer models selects particular subsets 
of four, two and one columnar dimer coverings as ground 
states: The '4-GS' model favors the four columnar states in 
X and y directions, e.g. states 3 — 6 in Fig.[T] and is described 
by the Hamiltonian 



7^4- 



GS 



V/) 



(2) 



The '2-GS' model favors columnar orderings only along one 
lattice direction, say the z direction, e.g. states 1 — 2 in Fig.[T] 
with Hamiltonian 



H2- 



GS 



^11 • 



(3) 



While these two models have somewhat different ground- 
state properties, they will turn out to exhibit rather similar 
physics. More specifically they share the same set of sym- 
metries, which we will discuss in detail in section [rV A2| The 
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last member of this first family of models is the '1-GS' model 
which singles out one of the six columnar orderings as sole 
ground state. We choose one of the columnar orderings in 
the z direction and define n|^|^° to be the number of plaque- 
ttes with parallel dimers on neighboring even or odd bonds, 
respectively. Then the Hamiltonian of the 1-GS model can be 
written as 



(4) 



A common feature of all four models introduced above is 
that they all undergo a direct thermal transition between the 
Coulomb gas phase at high temperature and a conventional 
long-range ordered state at low temperature as we will dis- 
cuss in Sec.|ll] Our parent Hamiltonian, the 6-GS model, has 
first been studied in Ref. 4, where based on an extensive nu- 
merical analysis the authors argue that this model undergoes a 
continuous thermal transition with the system spontaneously 
selecting one of the six columnar ordering patterns at the tran- 
sition out of the Coulomb phase. Our present numerical anal- 
ysis for the other three models indicates that the order of their 
respective transitions is (strongly) first order for the 2-GS and 
4-GS models, while we find strong evidence that the transition 
of the 1-GS model is again continuous. 

In Sec. [Ill] we provide further numerical evidence for the 
continuous nature of the 1-GS model by embedding this tran- 
sition into a line of continuous transitions for systems of iden- 
tical symmetry. The latter is achieved by studying a continu- 
ous interpolation of the 1-GS and 2-GS model in Sec. |III A| 
A similar approach interpolating the 4-GS to 6-GS model is 



given in Sec. Ill B 



To understand the different phase transitions found numer- 
ically, we develop candidate field theories in Sec.|lV] To this 
end, we first rewrite the dimer models in terms of a compact 
U(l) gauge theory. A subsequent duality transformation then 
allows to make the monopole excitations in this gauge theory 
explicit and describe the phase transitions as Higgs confine- 
ment transitions driven by the condensation of monopoles. 
Finally, the candidate field theories are derived in terms of 
two complex fields (a CP^ field) coupled to a U(l) gauge 
field. The symmetry-mandated Landau-Ginzburg actions for 



the various dimer models are given in Sec. IV B and impli- 
cations for the phase transitions discussed. In particular, this 
analysis suggests that the 1-GS model undergoes a continu- 
ous transition in the 3D inverted XY universality class, while 
the 2-GS model should undergo a first-order transition. While 
analytically inferring the nature of the transition for the 4-GS 
and 6-GS model is somewhat more delicate as discussed in 



Sec. IV B| we find an overall good agreement with the numer- 
ical results. 

In Sec. [V]we discuss a second family of dimer models with 
the common characteristic that there are two subsequent ther- 
mal transitions between the Coulomb phase and the dimer 
crystal. As a representative model we discuss the so-called 
'xy' model in some detail, which favors one columnar order- 
ing pattern along the x and y lattice directions, respectively. 
We find that the system first undergoes a continuous Higgs 
transition out of the Coulomb phase, which again is in the 3D 



inverted XY universality class, and subsequently a first-order 
"spin-flop" transition to the dimer crystal. 

We give some concluding remarks in Sec.|Vl]on the general 
interest of these classical dimer models in the search of non- 
LGW phase transitions. 



n. MONTE CARLO SIMULATIONS 
A. Overview 

We first summarize some characteristic numerical results 
for the thermal transitions obtained from extensive Monte 
Carlo simulations. The classical nature of our models not only 
commands to use an efficient stochastic algorithm to traverse 
the space of dimer coverings, but also allows for non-local 
update schemes such as the worm algorithm'"*'^. The latter 
performs an update by flipping a whole sequence of dimers 
when moving from one dimer covering to another one, which 
drastically reduces the problem of critical slowing down close 
to phase transitions. We used this algorithm on samples with 
N ^ I? sites up to 256^ - sizes which sometimes turn out to 
be necessary to ascertain the nature of a phase transition. 

The update scheme of the worm algorithm further allows to 
sample the behavior of two test monomers embedded in the 
dimer coverings. In particular, the update is performed by ini- 
tially breaking up an arbitrary dimer into a pair of monomers 
and then moving one monomer across the lattice by flipping 
dimers along a string or 'worm' until it can be recombined 
with the other monomer into a newly formed dimer. This con- 
struction can be used to reveal the confining properties of the 
low -temperature phases in our dimer models. To this end, 
we define the monomer 'confinement length' ^{T) as the 
(squared) average distance between the two test monomers, 
which we rescale by the expectation value {I? + 2) /4 for de- 
confined monomers moving freely on the lattice for a finite 
cube of even linear extent L (and periodic boundary condi- 
tions). 

We also measure thermodynamic al quantities such as the 
internal energy E{T), the specific heat C„(r) = ((i?^) - 
{E)'^)/T'^ as well as the stiffness p. The stiffness encodes 
fluctuations of dimer fluxes: p — y zi^l,) /'■^^^ where 

the flux is the algebraic number of dimers crossing a plane 
perpendicular to the unit vector Bq. Algebraic here means 
that, given a lattice direction, we count +1 for a dimer going 
from one sublattice to the other and —1 for the reverse sit- 
uation. Fluxes are conserved quantities (plane by plane) 
which vanish on average for symmetry reasons. 

In Fig.|2]we plot the specific heat per site C„ (T) /N, the en- 
ergy per site E{T) /N, and the monomer confinement length 

(r) for the four models introduced in the previous section. 

For all four models we find clear thermodynamic signatures 
for a direct transition between the high-temperature Coulomb 
phase (with deconfined monomers) to the dimer crystal at low 
temperatures (with confined monomers). The sharp, kink-like 
features in the energy E{T) and monomer confinement length 
^^(r) are indicative of a first-order transition for the 2-GS 
and 4-GS models. The first-order nature of these transitions 
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FIG. 2: (color online) Overview of numerical results for the dimer 
models with one (1-GS), two (2-GS), four (4-GS) and six (6-GS) 
columnar ground states: a) the specific heat per site Cv(T)/N, b) 
the energy per site E{T) /N, and c) the monomer confinement length 
^^(r) defined in the text. All four models undergo a direct transition 
with clear thermodynamic signatures between the high-temperature 
Coulomb phase (with deconfined monomers) to the dimer crystal at 
low temperature (with confined monomers). The sharp, kink-like 
features in the energy E{T) and monomer confinement length ^^(T) 
are indicative of a first-order transition for the 2-GS and 4-GS mod- 
els. The 1-GS and 6-GS model appear to undergo continuous transi- 
tions with smooth features. Data is shown for system size L = 48 
for the 1-GS and 6-GS models, and L = 32 for the 2-GS and 4-GS 
models. 
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FIG. 3: (color online) Specific heat per site Cv{T)/N for the 1- 
GS model as function of temperature T and different system sizes 
L. The inset shows a specific heat scan over a wide temperature 
range for a sample L = 48 system size. Below the peak around 
Tc ~ 2.276 ± 0.001 there is a shoulder which shows no variation 
with system size. 



is fully revealed in bimodal energy histograms in the vicinity 
of the transition temperature, which we will discuss in detail 
in section 



111 The 1-GS and 6-GS appear to undergo contin- 



uous transitions with smooth features in the energy E(T) and 
monomer confinement length ^■^ (T) and a (divergent) peak in 
the specific heat Cy{T). We have found no evidence of bi- 
modal energy histograms in the vicinity of these transitions, 
as also discussed in sectionlllll 

We will now turn to the individual dimer models and dis- 
cuss our numerical results in more detail in the following. 



B. The 1-GS model 

We will first concentrate on the 1-GS model which energet- 
ically favors a single columnar dimer ordering pattern shown 
in Fig. [T] Our numerical simulations for systems with up to 
256'^ dimers clearly suggest that this model undergoes a con- 
tinuous thermal transition between the Coulomb phase and the 
dimer crystal. 

The specific heat plotted in Fig. [5] exhibits a peak around 
the transition temperature of Tc ~ 2.276 ± 0.001 that appears 
to diverge very slowly with L. Below this peak there is a 
shoulder that does not show any variation with system size, 
(see inset of Fig. [3]l, and thus cannot be associated with any 
long distance or critical behavior. The latter is reminiscent 
of the 6-GS model"* which below the transition temperature 
exhibits an even more pronounced shoulder (for a comparison 
see also Fig.|2]i. 

A distinct feature of the Coulomb phase is that (test) 
monomers are deconfined. As a consequence, we expect 
the monomers to confine at the phase transition out of the 
Coulomb phase. This confinement transition can be tracked 
using the monomer confinement length £,'^{T) introduced 
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FIG. 4: (color online) Confinement length (T) for the 1-GS model 
measuring the (squared) average distance between two monomers. 
Data for different system sizes L are renormalized by the expecta- 
tion value (L'^ + 2)/4 for deconfined monomers. The distinct cross- 
ing point shown in the inset for temperatures in the vicinity of the 
transition temperature indicates a continuous transition. 
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FIG. 5: (color online) Stiffness p of the 1-GS model multiplied by 
L versus temperature T near the transition, for different system sizes 
L. The crossing point indicates a continuous transition. 



above. Plotting data for various systems sizes, as shown in 
Fig. [4] we observe a distinct crossing point at the transition 
temperature. This absence of finite-size effects at the tran- 
sition temperature indicates a universal value of the confine- 
ment length ^(Tc) at this transition, which we estimate to be 
^(Tc) ~ 0.923 ±0.001. This crossing point strongly indicates 
a continuous transition. 

Another indication of a continuous transition is that the 
distribution of dimer fluxes also becomes universal at the 
transition temperature. Indeed we observe a distinct crossing 
point for the stiffness p (multiplied by system size L) when 
plotting curves for different L in the vicinity of the transition 
out of the Coulomb phase, as shown in Fig.|5] The position of 
this crossing point coincides exactly with the transition tem- 



2.4 
2 
1.6 

Q. 

» 1.2 
S 0.8 
0.4 




♦ L = 256 




■o. 


L = 224 






= L=192 






L=160 






= L=128 
= L=112 
° L = 96 




L = 80 






« L = 64 




/ ; 


L = 48 






> L = 32 
: " L = 24 


y 


V = 0.6717 (3D XY) : 


- • L=16 




T =2.27618 ^ 






C 


1 


1 ' 1 


111 



2.4 
2 

1.6 
1.2 
0.8 
0.4 




, 1 /V 

L t 



00 0.9 



0.8 - 



0.7 



L= 128 
L = 96 
Z, = 64 
Z, = 48 
L = 32 
L = 24 
L= 16 



v = 0.6717 (3D XY) 

T =2.27618 



0.9 



- O.S 



0.7 



,1 /V 

L t 



FIG. 6: (color online) Data collapse for the stiffness p multiplied 
by L (top panel) and confinement length (lower panel) of the 1- 
GS model as a function of L^l^t, where t = (T ~ Tc)/Tc with 
Tc — 2.27618. The critical exponent = 0.6717 corresponds to the 
3D XY universality class. 



perature ~ 2.276 ± 0.001 estimated from the specific heat. 

Having established the continuous nature of the transition, 
we now turn to its universality class. Since this phase transi- 
tion occurs without any spontaneous symmetry breaking, we 
cannot rely on conventional techniques using an order param- 
eter to measure critical exponents. However, we can still con- 
sider thermodynamics, such as the behavior of the specific 
heat in the vicinity of the transition. As shown in Fig. |3] 
Cv{Tc)/N grows very slowly with system size at criticality, 
which would suggest a critical exponent a > 0, but very 
small. It is also quite possible that C„ (T^) /N actually con- 
verges to a finite value, but for system sizes that are currently 
out of reach of our numerical simulations. This would indi- 
cate a negative value for a < 0, also likely very small. This 
latter scenario is not unlikely considering the 3D XY model, 
which is known to have a small negative critical exponent 
a = —0.0151'^, but for which numerical simulation do 
not see a convergence of the specific heat. 

Thermodynamics being of little help to determine the uni- 
versality class, another possibility is to consider crossings and 
data collapse of adequate quantities, including the stiffness 
and the confinement length. Standard finite-size scaling argu- 
ments indicate that close to the transition point, the stiffness 
should scale as p = j^p{L^^'^ ■ t), where p is a universal func- 
tion, t = {T — Tc)/Tc the deviation from the critical tem- 
perature, and u the correlation length exponent. Performing 
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FIG. 7: (color online) Confinement length ^^{T) for the 2-GS 
model. Data for different system sizes L are renormalized by the ex- 
pectation value (L^ + 2) /4 for deconfined monomers. The absence 
of a distinct crossing point in the vicinity of the transition tempera- 
ture indicates a first-order transition. 



FIG. 8: (color online) Confinement length ^^(T) for the 4-GS 
model. Data for different system sizes L are renormalized by the ex- 
pectation value (i^ + 2) /4 for deconfined monomers. The absence 
of a distinct crossing point in the vicinity of the transition tempera- 
ture indicates a first-order transition. 



this analysis, we find a nice data collapse for the correlation 
length exponent v = 0.6717 of the 3D XY universality clas^ 
as shown in the top panel of Fig. |6] The same scaling form 
^ = i{L^^'^ ■ t) is also expected for the confinement length 
As shown in the lower panel of Fig. [6]we again find a data col- 
lapse for the same exponent v — 0.6717. Finally, we note that 
the system-size independent value ^(T^,) w 0.923 ± 0.001 is 
another characteristic of the universality class of the transition 
and in this case also points to the 3D XY universality class 



C. The 2-GS and 4-GS models 

We contrast our finding of a continuous transition in the 1- 
GS model with some numerical results for the 2-GS and 4-GS 
models which both undergo first-order transitions between the 
Coulomb phase and the dimer crystal phase. In Figs. |7]and[8] 
the confinement length of monomers is plotted. Similarly to 
the 1-GS model the transition out of the Coulomb phase is 
accompanied with a confinement of the monomers. At these 
first-order transitions we do not observe a distinct crossing 
point, see the insets of Figs. |7] and [8] 



III. CONTINUOUS INTERPOLATION BETWEEN DIMER 
MODELS 

One way to firmly establish the continuous nature of the 
phase transition in the 1-GS and 6-GS dimer models is to 
demonstrate that these transitions can be embedded into lines 
of continuous transitions. We will first concentrate on the 1- 
GS model and show that such a line of continuous transitions 
ending in a proposed, multicritical point can indeed be found. 
We will then discuss a similar idea for the 6-GS model, which 
however does not reveal such a line of continuous transitions. 



A. Interpolating the 1-GS and 2-GS models 

We have already observed that the 1-GS model which fa- 
vors a single columnar ordering pattern in one lattice direc- 
tion undergoes a continuous transition, while the 2-GS model 
which favors the two possible columnar ordering patterns in a 
given lattice direction undergoes a strong first-order transition. 
We can now ask how the nature of the phase transition changes 
as we continuously interpolate between these two models. To 
this end we continuously vary the weights for the columnar 
ordering patterns on the odd/even bonds in a given lattice di- 
rection. Formally, we introduce a coupling parameter A with 
< A < 1 on every other bond 



7^i_2-GS--^(Anf|+nf|) 

□ 



(5) 



For A = we recover the 1-GS model, while A = 1 corre- 
sponds to the 2-GS model. 

Our numerical simulations for various interpolation param- 
eters A are summarized in Figs.|9]and 10 which show the en- 
ergy per site E{T) /N and histograms of the energy per site in 
the vicinity of the transition temperature, respectively. Start- 
ing from the 2-GS model (A = 1) the sharp, kink-like feature 
in the energy accompanying the first-order transition quickly 
vanishes for interpolation parameters A < 1 as the two possi- 
ble columnar dimer orderings acquire different weights. The 
energy histograms in the vicinity of the transition tempera- 
tures turn from a bimodal distribution in the parameter regime 
1 > A ?i 0-97 into a single peak distribution for A < 0.97 



and system size L — 48, see Fig. 10 This strongly suggests 



that the first-order transition of the 2-GS model turns contin- 
uous for some intermediate A, which for larger system sizes 
might be closer to Ac ~ 0.95. On the other hand, this demon- 
strates that the continuous transition of the 1-GS model is in- 
deed part of a line of continuous transitions which extends 
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FIG. 9: (color online) Energies per site for a dimer model which FIG. 11: (color online) Energies per site for a dimer model which 
continuously interpolates between the 1-GS (A = 0) and 2-GS (A = continuously interpolates between the 4-GS (A = 0) and 6-GS (A = 
1) models. Data shown is for system size L = 48. 1) models. Data shown is for system size L = 48. 
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FIG. 10: (color online) Energy (per site) histograms for the dimer 
model interpolating between the 1-GS and 2-GS models in the vicin- 
ity of the transition temperature of the respective models. Data 
shown is for system size L = 48. 



over the range < A < Ac and likely ends in a multicritical 
point at Ac- Note that while A 7^ 1 introduces a staggering 
with respect to the two columnar ordering patterns, the sys- 
tem exhibits identical symmetries for all < A < 1, which 
will lead to a uniform theoretical description of the interpo- 
lated models for all < A < 1. However, since the strong 
first order transition of the 2-GS is expected to be stable to- 
wards a small perturbation (5A = 1 — A it is not surprising 
to see that the interpolated models exhibit (weak) first-order 
transitions in the regime 0.97 < A < 1, but quickly turn to 
uniform behavior for smaller A. 



B. Interpolating the 4-GS and 6-GS models 

We now turn to the 6-GS model, for which we can study 
a similar interpolation to the 4-GS model, which in contrast 
to the 6-GS model undergoes a first-order transition. Again 



we introduce an interpolation parameter A which now assigns 
different weights to the two columnar ordering patterns that 
establish the difference between the 4-GS and 6-GS model. 
Formally, we investigate the Hamiltonian 



n 



4-6-GS 



□ 



An||) , 



(6) 



with < A < 1. For A = this is the 4-GS model, while 
A = 1 now corresponds to the 6-GS model. 

Our numerical results for various interpolation parameters 
A are summarized in Figs. [TT| \V2\ and 13 respectively. We 
find that a sharp, kink-like feature around the transition tem- 
perature in the energy E{T) persists for almost all interpola- 
tion parameters < A < 1 (see Fig. 111. Energy histograms 



in Fig. 12 for the respective transition temperatures show bi- 
modal distributions for the same parameter range. Pushing the 
limit of our calculations we can establish a two-peak structure 
up to A — 0.98 for system size L — %Q (see the lower panel 
in Fig. 
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If we systematically trace the distance between 
the two peaks in these bimodal energy distributions, as shown 
in Fig. [13] the emerging trend clearly suggests that the line 
of first-order transitions persists all the way up to A < 1, but 
the transition turns weaker along this line with the data sug- 
gesting that there is no bimodal distribution for A = 1 and the 
transition becomes continuous for the 6-GS model. 



IV. CANDIDATE FIELD THEORIES 

To understand the different nature of the phase transitions 
discussed in Sec.|ll] we will now develop a family of candidate 
field theories describing these transitions. To do this, we take 
advantage of a pair of mappings: from the dimer model to 
compact lattice quantum electrodynamics (QED), and thence 
to a dual monopole formulation. The continuum limit of the 
latter leads directly to the desired field theories. 

The desired mappings have been discussed in some detail in 
Refs |19l3l and in fact can be performed not only for the clas- 
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FIG. 12: (color online) Energy (per site) histograms for the dimer 
model interpolating between the 4-GS to 6-GS models in the vicinity 
of the transition temperature of the respective models. Data shown 
are for system size L = 48 in the upper panel and L = 96 in the 
lower panel. 



sical dimer model discussed here but also its quantum gener- 
alization. For clarity, we will present in this section a brief, 
self-contained summary of the mapping in the classical case 
for the cubic lattice, applicable to the models of the present 
paper Performing a detailed symmetry analysis, we then de- 
rive the symmetry allowed Ginzburg-Landau actions for the 
various models and discuss implications for thek respective 
phase transitions. 



A. Mappings 

To proceed, we define a bond variable which counts the 
number of dimers on a given bond: 

f 1 , if the bond (a, b) is occupied by a dimer 

\0 otherwise. ^ ' 

The close-packed dimer constraint which requires that every 
site in the cubic lattice be part of exactly one dimer can then 
be expressed as 



b 



(8) 



FIG. 13: (color online) Distance between the two peaks in the bi- 
modal energy histograms shown in Fig. [12] The lines are guides to 
the eye. Data for system sizes L — 48 and L = 96 are shown. 



where the sum is over sites b which are nearest neighbors of 
a. A monomer excitation in the dimer model which breaks the 
close-packing constraint, e.g. an unpaired site a on the cubic 
lattice, is then indicated by Na = 0. 



1. compact QED 

We may directly pass to QED variables as follows. We in- 
troduce an electric field variable Eab, which is a directed vari- 
able, according to 



Eab — ^aflab 



(9) 



where Eab is integer-valued (in particular Eab = 0, ±1) and 
we have introduced a 'background charge' with a fixed dis- 
tribution of alternating charges on the two sublattices 



, a e A sublattice, 
— 1 , a e B sublattice. 



(10) 



In the QED formulation, the local constraint ([8]l maps di- 
rectly to a lattice version of the Gauss law. 



div£; 



(11) 



which also explains the notion of the background charge and 
where we have used the lattice divergence div E = J2b ^ab- 
Expressed in the QED variables, the Hamiltonian becomes 

Wqed 



□ 



2 ^ 

(a,b) 



curl E, 2 



const. 



curl £,-2/ 

(12) 



where the first term is a constant in the physical space in which 
hab = 0,1. We include it, however, in order that we may al- 
low the electric variable Eab to fluctuate over all integers; by 
taking the large U limit, the physical dimer states, which min- 
imize this term, are selected. It is expected that the universal 
properties of TYqed are identical for infinite and finite U . 



Note that the Hamiltonian ( [T2| i is rather similar to the stan- 
dard formulation of compact QED. The main difference is 
the absence of any magnetic field terms B^, which reflects 
the classical nature of the dimer model under consideration, 
a shifting of the term by an alternating 'background field' 
£0/2 and the energetic preference for curl£? = ±2. Despite 
these differences, Hamiltonian ( [T2j i does share all the same 
internal symmetries as the more conventional QED form. It 
is therefore expected to share the same properties in regimes 
where universality is mandated. 



2. Duality, Monopole Formulation and Symmetry Transformation 

In contrast to a conventional, non-compact QED formu- 
lation a compact QED like the one introduced in the previ- 
ous section does not prohibit magnetic monopoles. These 
monopoles are 'conjugate' to the gauge charges in the QED 
description, which in terms of the original dimer models cor- 
respond to monomer excitations. The defining characteristic 
of the Coulomb phase is that the gauge charges are decon- 
fined, while the magnetic monopoles are well-defined, gapped 
quasiparticles. In the dimer crystal phase, on the other hand, 
the electric charges are confined and the electric field is static. 
This implies that the conjugate magnetic field is strongly fluc- 
tuating and the magnetic monopoles are no longer good ex- 
citations. We can thus describe the phase transition out of 
the Coulomb phase as the (Bose) condensation of magnetic 
monopoles accompanied by a simultaneous confinement tran- 
sition of the electric charges. 

Our goal now is to establish an analytic description of this 
phase transition in terms of the magnetic monopoles. To this 
end, we will first introduce a duality transformation to make 
the monopole excitations explicit in the Coulomb phase. It is 
well-known that the electric and magnetic fields in Maxwell's 
equations are dual in the absence of charges and currents. 
However, while there are no currents in our system, there 
is a non-vanishing charge distribution, as given by the non- 
vanishing electric field divergence in Eq. dTTJ. We therefore 
introduce a 'background electric field' e '■%which compen- 
sates for these background charges by satisfying 



div e 



(13) 



and making i? — e divergence free. Note that while i? is a fluc- 
tuating variable, the background field e is static. It is conve- 
nient to choose e to be integer valued, for what follows. Then 
a simple choice to satisfy ( 13 1 is to place Ca^a+x = (1 + ^a) /2 



on the link emanating in the x-direction from site a. 

This allows us to write down a duahty transformation of the 
form 



E. 



ah 



curl d + (iab 



(14) 



where we have introduced a dual electric vector potential a. 
In a quantum theory, a would generate the conjugate magnetic 
field. Due to the integer constraint on Eab, we must also take 
a to be integer valued (this reflects our integer choice of Cab)- 




FIG. 14: (color online). Cubic lattice and its dual. The dual cubic 
lattice sites sit at the centers of each cube of the cubic lattice. 



In terms of these dual variables, the QED Hamiltonian ( fT2| 

becomes 



"^dual-QED 



U 



curl Q? — e 



(ab) 



(15) 



where the squares '□' now denote the plaque ttes of the dual 
cubic lattice, see Fig. 14 and the last term, i^[V^Q;] represents 



the transcription of the last term in Eq. ^V2\ in terms of a. 
We will not need its explicit form here. It is only important 
to note that it contains second order (lattice) derivatives of 
the vector potential a. Upon coarse-graining, such terms are 
irrelevant in the continuum limit. What is important is that, 
in the process of integrating out short scale fluctuations, they 
will generate relevant terms of all possible types dictated by 
symmetry. In this way, the physics of the dimer interactions, 
which is reflected in the F function, enters the low energy 
continuum field theory description. 

We now proceed to develop the continuum limit, following 
a sequence of standard manipulations^^. In doing so, we will 
neglect the F term, on the grounds discussed above, keeping 
in mind that in the final continuum theory, we must restore 
all possible symmetry-allowed interactions that may be gen- 
erated from it. We first soften the integer constraint on the 
vector potential a, replacing the constraint by a term which 
favors integer values. This approximation does not change 
the nature of the monopole condensing phase transition. We 
rewrite the Hamiltonian ( [T5| ) as 



r,r' 



COS [27 



(16) 



where large w recovers the integer a constraint. With this 
rewriting, we may regard a as a real-valued variable. In the 
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first term we have also dropped the 6^/2 term found inside 
the parenthesis of the first term in Eq.([T5]l. This is possible 
because this term, regarded as a vector field, is purely longitu- 
dinal (i.e. curl-free), and hence, actually decoupled from the 
curl a factor. We will soon extract a further longitudinal piece 
from e. 

To proceed, we first introduce explicit monopole phase 
variables by making the gauge transformation arr' curr' — 
{6r — 9r>)/2TT. One obtains 

7i = — (curia — e)"^ — w cos[9r — 0^' — 2TTary] ■ 

□ r,r' 



longitudinal parts. 



e = eL + er, 



Next we take the continuum limit, following Sec. VI. B of 
Ref I9j The hopping Hamiltonian in Eq. ( |22] i has two min- 
ima, which we will call ^'i(r) and 'i'2{r) here. A solution of 
the tight-binding Hamiltonian then becomes a linear combi- 
nation 



^(r) = (?!>i*i(r) + 02*2(r), 



(23) 



such that div eV = and curl = 0. Note that, because 
of the curl in Eq. (Vf) , only eV couples to a. Taking the di- 
vergence of Eq. (|18|l, we see that divet = dive = ea- A 
choice for satisfying this condition and which is curl-free 
is simply 

(eL)a,a+M - ea/6. (19) 

From this, we of course can find er by solving Eq. ( fTSj ): 



(20) 



Inserting Eq. ( fTSj ) into Eq. ( [T7] i, and dropping the decoupled 
and constant e l part, we find 

7i = — (curia — eV)^ — w cos[0,. — 9,.' ~ 27rar,r'] • 

□ r,r' 

(21) 

At this point, we have obtained a lattice Ginzburg-Landau 
theory, in which eV appears as the (average) dual flux (expe- 
rienced by the monopoles) through the dual plaquette pierced 
by this vector, expressed in units of the flux quantum. As 
usual, only the fractional part of the flux has physical signifi- 
cance. From Eq. ( |20l ), we readily see that this fractional part 
is uniformly 1/6 of a full flux quantum piercing the dual pla- 
quettes emanating from one sublattice of the direct lattice and 
ending in the other. This can be seen as an array of alternating 
dual monopole fluxes, representing the original alternating 
background charges in the QED theory. 

Precisely this problem, of monopoles moving in this back- 
ground flux pattern on the cubic lattice, was studied by 
Motrunich and Senthil, in Ref. [19j We can adapt their results 
directly. We define a soft-spin monopole field V'r ~ e*^'', 
and neglect at first the fluctuations in the gauge field, replac- 
ing arr' by a static gauge configuration cfrr' representing the 
background flux. The soft-spin Hamiltonian is 



monopoles 



h.c] . (22) 



where we treat 0i and 02 as slowly varying fields. 

We can now ask how these solutions transform under the 
symmetry operations of the original cubic lattice. For the spe- 
cific gauge choice of Ref.[19|these were reported to be 



(24) 



where the two 90° rotations RiL ji^y and RK^ri.xz along the 
z and y lattice directions are around the sites on which the 
monopoles reside - the center of the cubes of the original cu- 



(17) 


r : 


01- 


-0^ 


02 - 


--0^ 


into transverse and 


Ty: 


01- 


-01, 


02 - 


- 02 


< 


: 


01- 




02 - 


-0*1 




R^.Rxy ■ 


01- 




02 - 


e'"^/40; 


(18) 


R^,B.xz ■ 


01 - 


-75(0*1 + 02), 


02 - 





bic lattice / the sites of the dual lattice, see Fig. 14 

We have now established the symmetry transformation 
properties for the slowly varying fields in the solution of the 
gauge mean-field Hamiltonian (|22]). As we will describe in 
the next section this allows us to make an explicit connec- 
tion of these solutions to the individual members in our fam- 
ily of dimer models. In particular, we directly show how the 
symmetries of the microscopic dimer interaction, implicitly 
contained in the appropriate F term in Eq. ( fTSj l, re-enters the 
continuum field theory. 



B. Effective Ginzburg-Landau actions and pliase transitions 

We will now turn to the individual members in our family 
of dimer models and derive an effective description in terms 
of a Ginzburg-Landau action that respects the symmetries of 
the various models. This action is typically given in terms of 
the two slowly varying complex fields 0i and 02 coupled to 
the dual U(l) gauge field a^. We then analyze the derived 
actions and discuss the nature of the phase transitions in these 
field theories. 

Let us first establish some notations and introduce a three- 
component vector iV, which will serve as an order parameter 
indicating which dimer ordering is chosen as ground state 



^^(01,02) = 0l'?Q/30fi 



(25) 



where a are the three Pauli matrices. The six columnar ground 
states of our dimer models depicted in Fig. [T|then correspond 
to iV pointing along positive or negative x, y, z directions, re- 
spectively. Finally, let = (0i, 02) be a two-component vec- 
tor combining the two complex fields 0i and 02. 

We can now write the effective Ginzburg-Landau action as 



5eff = f3 d^r[ (V - id)0 + U{N) + C{a) , (26) 
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dimer model 


symmetries 


6-GS 


Tx, Ty, Tz, R^^Rxy, R^,Rxz 


4-GS 


Tx , Ty , Tz, R^Rxy 


2-GS 


Tx , Ty , Tz , ,Rxy 


1-GS 


Tx , Ty , ^Rxy 


xy 


Tz , R'^^Rxy 


xxy 


Tx , Tz 


xyz 


R^,Rxy, R^,Rxz 


xyzz 


Tz, R^,Rxy 


xxyyz 


Tx, Ty, R:^^Rxy 



TABLE I: Symmetries of the various dimer models. 



where the first | . . . p-term is a minimal coupling of the two 
complex fields to the dual U(l) gauge field and £(q?) is the 
usual Maxwell's term for the gauge field. The potential U (N) 
is determined by the underlying symmetries of the various 
dimer models, which are summarized in Table [l] This poten- 
tial therefore varies for the individual models as discussed in 
more detail the following. 

Note that the action ( p6] l does not contain any time- 
derivatives. The reason is that in the presence of such time- 
derivative terms and periodic boundary conditions in imagi- 
nary time all modes with non-zero Matsubara frequencies are 
more massive than the zero frequency mode of interest here 
and can be integrated out. 

1. The 1-GS model 

The ground state of the 1-GS dimer model is a single 
columnar ordering pattern, which we choose to be oriented 
along the z direction. The effective Ginzburg-Landau action 
for this model is thus required to be invariant under the trans- 
formations Tx, Ty and RiL j^xy only. In particular, the poten- 
tial U {N) in the action (|26| has the general form 

(27) 

where we have only included terms up to quadratic order in 
the complex fields and introduced two coupling constants U2 
and V2- Note that the potential U (N) introduces two inequiv- 
alent mass terms for the two complex fields. As we reduce 
temperature it will be the complex field with smaller mass 
that will condense thereby leading to a condensation of the 
monopoles (while the other field still has vanishing expeca- 
tion value). The corresponding phase transition can thus be 
described by a field theory with just one complex field cou- 
pled to a U(l) gauge field which is known to be a continu- 
ous transition in the inverted 3D XY universality class^'. At 
this Higgs transition the system spontaneously breaks the U( 1 ) 
gauge symmetry, but does not break any lattice symmetries. 
As the Coulomb phase breaks down the charge excitations, 
e.g. the monomers in the language of the dimer model, con- 
fine. 

We thus have clear analytical and numerical evidence for a 
continuous transition between the Coulomb phase and a long- 



range ordered dimer crystal in this dimer model, which cannot 
be explained by the standard LGW paradigm. 

2. The 2-GS and 4-GS models 

We now turn to the 2-GS and 4-GS models which we have 
seen to exhibit direct first-order thermal transitions. Although 
the two models have complementary ground-state manifolds, 
they are invariant under the exact same lattice transformations. 
As given in Table [ijthese are the three lattice translations T^, 
Ty, and as well as rotation around the z-axis, RiL^nxy The 
symmetry allowed potential terms up to quartic order in the 
Ginzburg-Landau action (|26| are thus given by 

U{N)^U2\N\+Ui\N\'^ + vaNI . (28) 

For the 4-GS model U4 > 0, so the term modulates 
by the constraint \(f)i\ — |(/)2| and prefers the order parameter 
N to point in the xy plane. Noteworthily, this is still a con- 
tinuously connected manifold with an internal U(l) symmetry 
for the order parameter Note that the system does not need 
to break any lattice symmetries to satisfy this constraint (in 
contrast to the 6-GS model which we will discuss in the next 
section). At the transition, when the complex fields (f) con- 
dense, the order parameter N becomes non-zero and points 
along one of the four lattice directions in the xy plane. Note 
that at this Higgs transition the system not only spontaneously 
breaks the U(l) gauge symmetry, but simultaneously also the 
U(l) order parameter symmetry as well as the four-fold lattice 
symmetry. 

The action ( |28| ) exactly corresponds to the one studied 
for an easy-plane quantum antiferromagnet in the context of 
deconfined quantum criticalitjB^EU While analytical inves- 
tigations of this action have suggested a continuous phase 
transition^^, extensive numerical results have pointed to a 
weak first-order transitio n^^P^ , which is also what we find in 
our present numerical analysis. 

For the 2-GS model the order parameter N wants to point 
along the z direction, e.g. becomes maximal, which im- 
plies (contrary to the 4-GS model) that 114 < 0. This leaves 
the system with a disconnected manifold (of two points) for 
fixed magnitude \N\ either prefering |0i| = or |02| = and 
resulting in a Z2 symmetry for the order parameter. 

In contrast to the 1-GS model, this theory cannot be reduced 
to a field theory with just one complex field without breaking 
the lattice symmetry Tz- One possibility now is to have two 
subsequent transitions where we first break the lattice symme- 
try at a higher temperature and subsequently observe a Higgs 
transition at a lower temperature (with an exotic intermediate 
phase of coexisting Coulomb and dimer crystal correlations). 
In terms of the complex fields </) the system would sponta- 
neously select one of the two possibilities 7^ 0, |(/i2| = 
or 101 1 = 0, 1^2 1 7^ at the first transition. At the second tran- 
sition the non-vanishing 0-field would require a fixed phase in 
a Higgs transition. 

Another possibility is to have one direct transition. How- 
ever, it is hard to imagine a field theory giving rise to a contin- 
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uous transition where the spontaneous breaking of the discrete 
Zi order parameter symmetry occurs simultaneously with the 
Higgs transition breaking the U(l) gauge theory. Thus, we 
conclude that a direct transition is likely first-order. It appears 
to be the latter scenario that we observe in our numerics. 

Since we have only numerical (and not analytical) evidence 
for a first-order transition in the 4-GS model, but have some 
analytical evidence for a first-order transition in the 2-GS 
model, we probably expect the latter to be the stronger first- 
order transition. This is what we observe in the numerical 
simulations of SecHIl 

3. The 6-GS model 

Finally, we turn to the 6-GS model which respects all the 
cubic lattice symmetries. In writing down a symmetry al- 
lowed potential for the action ( p6] l we consider the simplest 
invariants of the form 

U{N) = V{\N\) + vsh{N), (29) 

where /8(^) = N^N^ + N^N^ + N^N^ is an eighth order 

term in the complex fields (j) and we adopted a notation similar 
to the one of Ref.[T9] Again we can expand the potential term 
T/(|7V|) =U2\N\+U4\N\'^ + .... Omitting the 8th order term 
in ( |29] l gives an SU(2) invariant action, which has attracted 
some interest due to recent proposals of SU(2) invariant de- 
confined quantum critical points, as suggested in the J — Q 
quantum model^^ . 

Following the line of arguments in Ref. [19] the confining 
Higgs transition (which we observe) occurs for U2 < and 
|7V| simultaneously acquires a finite magnitude, e.g. the U(l) 
gauge symmetry and the lattice symmetry are broken at the 
same transition. As argued in Ref. one of the six colum- 
nar ground states is selected by the eighth order term ug/g (N) 
with vs > 0. Since the Higgs transition of the action (|26]l 
is suggested to be continuous without this 8th order term and 
this term is likely irrelevant (due to its high order), we con- 
clude that the action (|29| allows for a continuous transition. 
This seems to be in agreement with the numerical evidence 
of Ref. 4 and our present numerical analysis interpolating be- 
tween the 4-GS and 6-GS models. 

Direct numerical simulation of the SU(2) invariant action 
( p9l l without the 8th order term has provided controversial re- 
sults with some evidence for a continuous transition^^, while 
another recent analysis favors a weak first-order transition^^. 
Including the 8th order term in the action ( |29] l a direct numeri- 
cal simulation of a similar action to p9l ) reported a continuous 
transition with exponents close but apparently different from 
those of the 3D XY modeP. 



4. The interpolated models 

Finally, we briefly turn to the models interpolating between 
models exhibiting continuous and first-order transitions as dis- 
cussed in section Bill 



Interpolating between the 1-GS and 2-GS model the sys- 
tem exhibits for all interpolation parameters < A < 1 
the same lattice symmetries and is therefore described by the 
same Ginzburg-Landau action with potential ( p7j l as the 1-GS 
model. This symmetry analysis suggests that for all A < 1 
the two complex fields (pi, (p2 acquire different masses and 
we can describe the action in terms of a single complex field 
coupled to a U(l) gauge field. On the other hand, we expect 
the strong first-order transition of the 2-GS model (A = 1) to 
be stable towards small perturbations and therefore to extend 
over a finite region A < 1. As a consequence, there should 
be a multicritical point Ac where the line of continuous transi- 
tions for A < Ac meets the first-order line for A > Ac. 

Interpolating between the 4-GS and 6-GS model the system 
exhibits for all interpolation parameters < A < 1 the same 
lattice symmetries, while the symmetries change for the end- 
point A = 1 which corresponds to the 6-GS model, see also 
Table |l] Again this symmetry analysis suggests that all inter- 
polated models with < A < 1 are described by the same 
Ginzburg-Landau action with potential ( |28] l. This seems to 
be in agreement with our numerical results suggesting that all 
models with < A < 1 exhibit first-order transitions. 



V. AN INTERMEDIATE PARAMAGNET 

Finally, we turn to a second family of dimer models that 
also energetically favor specific subsets of the six columnar 
ordering patterns illustrated in Fig. [T] The distinct feature of 
this second family of dimer models is that they harbor two 
consecutive thermal phase transitions. The high-temperature 
phase transition out of the Coulomb phase is again driven by 
the condensation of monopoles with confining monomer exci- 
tations. However, this phase transition is into a paramagnetic 
phase without dimer crystalline order which only forms at the 
low-temperature transition. Thus, we are left with an unusual 
sequence of phases in these models with the paramagnet re- 
siding at intermediate temperature scales. 



A. A second family of dimer models 

Our second family of dimer models explores other com- 
binations of the six columnar ordering patterns in Fig. [T] as 
ground states. The common characteristic in selecting the ad- 
missible ground states is that for at least one lattice direction 
we choose only one of the two possible columnar orderings 
and there is more than one ground state. If we name the mod- 
els by the lattice directions for which ground states are chosen, 
these are the 'xy', 'xyz', 'xxy', 'xyzz' and 'xxyyz' models. 
In this nomenclature the models in our first family of models 
would be named 'z','zz','xxyy' and 'xxyyzz' for the 1-GS, 
2-GS, 4-GS and 6-GS model, respectively. 

We will not discuss all possible models in this section, but 
concentrate on the 'xy' model with Hamiltonian 

□ 
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FIG. 15: (color online) Overview of the xy-model that favors one 
columnar ordering pattern along the x and y lattice directions: a) 
specific heat per site Cv{T)/N, b) energy per site E(T)/N, and c) 
confinement length This model undergoes an unusual sequence 
of phases with a paramagnetic phase (shaded area) at intermediate 
temperature scales. At the high- temperature transition out of the 
Coulomb phase (denoted by an arrow) the monomers confine, while 
the dimer crystal forms only at the low-temperature phase transition. 
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FIG. 16: (color online) Crossing point of the confinement length 
^^(r) for different system sizes at the high-temperature phase tran- 
sition out of the Coulomb phase for the xy model. 



transition out of the Coulomb phase the monomer confine- 
ment length drops again indicating that this transition is due to 
monopole condensation. A finite-size scaling analysis reveals 
a distinct crossing point (see Fig.[T6]l indicating a continuous 
transition into the intermediate temperature paramagnet. We 
will argue that this transition is again described by the inverted 
3D XY universality class. 

We first notice that that the temperature of the Coulomb 
transition in the xy model {Tc ~ 2.247) turns out to be close 
to the one found for the 1-GS model (where Tc — 2.276). 
Another indicator that the Coulomb transitions in these two 
models are closely related and probably of the same universal- 
ity class is that the universal value of the confinement length 
at the crossing point is l{Tc) « 0.920 ± 0.001, which is 
rather close to the one found for the 1-GS model (^(Tc) « 
0.923 ± 0.001). As a final argument, we find an excellent 
data collapse of the confining length measured in the vicinity 
of this transition for different system sizes when rescaling the 
data with the correlation length exponent v = 0.6717 of the 



3D XY universality class, see Fig. 17 

At the low-temperature transition the system spontaneously 
selects one of the two possible columnar ordering patterns and 
we observe a sudden increase of the number of plaquettes with 
parallel dimers, e.g. or n'^^, respectively (see Fig. fTSi 



The sharp jump of this order parameter indicates a likely oc- 
curence of a first-order phase transition. The first-order nature 
of the low-temperature phase transition is indeed confirmed 
by the bimodal structure of energy histograms close to this 
transition point (see Fig 
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where we have chosen the columnar dimer orderings on the 
even bonds in the x and y lattice directions as ground states. 

We summarize our numerical results for this model in 
Fig. 
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The two consecutive thermal transitions both carry 
distinct thermodynamic signatures with a double peak struc- 
ture emerging in the specific heat. At the high-temperature 



B. Theoretical Analysis 

We can discuss the nature of the phase transitions for 
this second type of dimer models by again analyzing the 
symmetry-allowed effective Ginzburg-Landau actions in full 



analogy to the discussion in section IV B for the first family 
of dimer models. 
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FIG. 17: (color online) Data collapse for the confinement length ^ 
as a function of L^^^t, where t = {T - Tc) /Tc with Tc = 2.2475. 
The critical exponent v = 0.6717 corresponds to the 3D XY univer- 
sality class. 
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FIG. 18: (color online) Bond occupation along the x, y and z lat- 
tice directions versus temperature for the xy model. At low temper- 
atures data is shown for a system that spontaneously orders along 
the y direction (other systems spontaneously order along the x direc- 
tion). The low temperature transition between the paramagnet and 
the dimer crystal is accompanied by a sudden increase of the bond oc- 
cupation along the preferred lattice direction. The high-temperature 
transition between the Coulomb gas and the paramagnet where the 
monomers confine is indicated by the arrow. 



With the lattice symmetries for the individual members of 
this second family of dimer models given in Table |lj we find 
the following potentials for the Ginzburg-Landau action in 
Eq. (|26]) 



UiN) 



U2\N\ 
U2\N\ 
U2\N\ 
U2\N\ 
U2\N\ 



Ny), 



■V2Ny, 

-V2{N^ + Ny + N^), 

-V2{N^ 
-V2N,, 



Ny), 



xy model, 
xxy model, 
xyz model, (31) 
xyzz model, 
xxyyz model . 



For all models we can diagonalize the quadratic part in 
these potentials by performing a SU(2) rotation in the {(f>i , (1)2) 
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FIG. 19: (color online) Energy histogram in the vicinity of the low- 
temperature phase transition, which we argue to be a thermal analog 
of a spin-flop transition (see text). Our numerics indicate a weak first- 
order transition by a double-peak structure in the energy histogram, 
which emerges for larger system sizes. Here system size is i = 96. 



space such that the potentials ( |3T| take an identical form as 
given in Eq. ( p7j i for the 1-GS model. As a consequence, 
we expect the high-temperature transition out of the Coulomb 
phase in all these models to be described by the same Higgs 
mechanism we identified for the 1-GS model resulting in 
a continuous transition in the inverted 3D XY universality 
class. Our numerics give supporting evidence for a continu- 
ous transition in this universality class as discussed above. 

The main distinction between the models in our second 
family of models and those in the first set of models is that 
here we can break an additional lattice symmetry which ap- 
parently gives rise to the second transition into the dimer 
crystal phase at lower temperatures. We argue that this sec- 
ond low-temperature phase transition is generically a first- 
order transition analog to a spin-flop transition. To see this 
analogy consider the xy model with the potential U{N) = 
U2\N\+ V2{Nx + Ny) in the Ginzburg-Landau action. Below 
the confinement transition the order parameter N has a non- 
zero expectation value (since the monopoles are condensed), 
but points half-way between the x and y directions, thereby 
minimizing the second term in the potential. This is also ev- 
ident in our numerical simulations as shown in Fig. [18] At 
very low temperatures, however, we know that the system 
must (because there will be no dimer fluctuations) sponta- 
neously order along one of the two lattice directions, thus 
breaking the symmetry between x and y directions. There- 
fore the spin N must reorient away from the (110) axis to 
the (100) or (010) axis. To describe this, we require addi- 
tional higher-order terms in the potential U{N), of the form 
NxNy, N^Ny, etc. At low temperatures, since the magnitude 
of the spin becomes large (|iV| = 1 as there are no fluctua- 
tions), such terms are no longer negligible. On lowering the 
temperature and increasing these higher order terms, we ex- 
pect that the minimum directions of this energy function may 
abruptly switch to their low temperature values. This is in- 
deed the most commonly occurring situation in spin systems, 
in which such a first-order reorientation is known as a "spin 
flop" transition. This expectation, arrived at above from "an- 
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alytical" field theory considerations, is indeed verified in the 
numerics (see Fig.[T9]l. 

VI. DISCUSSION 

Recent years have seen an extensive search for continu- 
ous phase transitions beyond the LGW paradigm, which were 
originally suggested to occur in certain quantum modelJ^SSI 
In this manuscript, we have demonstrated that this exotic 
physics can manifest itself also in various classical models. 
This, of course, is not much of a surprise since the universal- 
ity of continuous phase transitions mandates that they occur in 
a large variety of models, including classical ones. Neverthe- 
less, it is amusing to note that such unconventional phase tran- 
sitions and the sophisticated ordering mechanisms associated 
with them can actually be found in simple variations of one of 
the golden models of statistical mechanics, namely the dimer 
model. The key ingredient giving rise to this exotic physics is 
a constraint which enforces close-packed coverings of hard- 
core dimers. In a way, this readily builds into the classical 
model a certain level of frustration which is often invoked to 
be a key ingredient for quantum models to exhibit non-LGW 
criticality. 

Besides the important step to directly establish the occur- 
rence of non-LGW transitions in these dimer models, we 
view several advantages arising from their classical nature: 
(i) Classical models are notoriously simpler to analyze, both 
theoretically and numerically, than quantum models. They 
are accessible to Monte Carlo approaches, thus allowing to 
study critical phenomena through the direct simulation of 



large systems. For the specific dimer models at hand, the 
existence of a highly-efficient Monte Carlo worm algorithm 
is also very attractive, (ii) These models ease the identifi- 
cation of the necessary ingredients that are needed to cause 
non-LGW physics in a lattice model (such as lattice and/or 
continuous symmetries). This further opens the possibility 
of 'reverse-engineering' or 'rolling back the path integral' 
to obtain two-dimensional quantum models that exhibit the 
same non-LGW criticality as their three-dimensional classi- 
cal counterparts. Such a classical-to-quantum mapping was 
recently used in Ref 31 (iii) The stability of certain critical 
behavior can be easily explored in variations of these classi- 
cal models, e.g. through the inclusion of perturbations or by 
extrapolating terms (as performed in the current study). For 
instance, a yet-to-be-explored possibility is to include terms 
that frustrate the columnar ordering. A similar situation in a 
quantum model would generically come with a sign problem 
in Quantum Monte Carlo simulations, putting serious limita- 
tions to any numerical study, (iv) Finally, the sheer simplicity 
of these models might indicate that non-LGW transitions are 
not that exotic after all. 



VII. ACKNOWLEDGMENTS 

Our numerical work used some of the ALPS librariePEl 
see also http://alps.comp-phys.org This work was supported 
by the DOE through Basic Energy Sciences grant DE-FG02- 
08ER46524. LB's research facilities at the KITP were sup- 
ported by the National Science Foundation grant NSF PHY- 
0551164. 



' D. Bergman, R. Shindou, G. Fiete, and Leon Balents, Phys. Rev. 

Lett. 96, 097207 (2006). 
^ D. Bergman, J. Alicea, E. Gull, S. Trebst, and L. Balents, Nature 

Phys. 3, 487 (2007). 
^ D. Bergman, G. Fiete, and L. Balents, Phys. Rev. B 73, 134402 

(2006). 

F. Alet, G. Misguich, V. Pasquier, R. Moessner, and J. Jacobsen, 
Phys. Rev. Lett. 97, 030403 (2006). 
^ T. S. Pickles, T. E. Saunders, and J. T. Chalker, EPL 84, 36002 
(2008). 

C. Castelnovo, C. Chamon, C. Mudry, and P. Pujol, Phys. Rev. B 

73, 144411 (2006). 
^ D. A. Huse, W. Krauth, R. Moessner, and S. L. Sondhi, Phys. Rev. 

Lett. 91, 167004 (2003). 
^ M. Hermele, M. P. A. Fisher, and L. Balents, Phys. Rev. B 69, 

064404 (2004). 

' S. V. Isakov, K. Gregor, R. Moessner, S. L. Sondhi, Phys. Rev. 
Lett. 93, 167204 (2004); C. L. Henley, Phys. Rev. B 71, 014424 
(2005). 

F H. Stillinger and M. A. Cotter, J. Chem. Phys. 58, 2532 (1973). 
" J. Villain and J. Schneider, in Physics and Chemistry of ice edited 
by E. Whalley, S. J. Jones, and L. W. Gold (Royal Society of 
Canada, Ottawa, 1973); see also J. Villain, Solid State Commun. 
10, 967 (1972). 

R. W. Youngblood and J. D. Axe, Phys. Rev. B 23, 232 (1981). 



" G. Misguich, V. Pasquier, and F Alet, Phys. Rev. B 78, 100402(R) 
(2008). 

A. W. Sandvik and R. Moessner, Phys. Rev. B 73, 144504 (2006). 
F. Alet, Y. Ikhlef, J.L. Jacobsen, G. Misguich, and V. Pasquier, 
Phys. Rev. E 74, 041 124 (2006) 

M. Campostrini, M. Hasenbusch, A. Pelissetto, and E. Vicari, 
Phys. Rev. B 74, 144506 (2006). 

" A Cucchieri et al, J. Phys. A: Math. Gen. 35, 6517 (2002). 

A similar value ^^(T^) = 0.9204(10) is found for the 3D 
isotropic link-current model (see for instance F. Alet and E. 
S0rensen, Phys. Rev E 67, 015701 (2003)) which has a phase 
transition in the 3D XY universality class (F. Alet, unpublished). 
The comparison with the current study needs a further investiga- 
tion of the role of anisotropy in the link-current model (as the 
dimer models are themselves anisotropic). 

" O. Motrunich and T. Senthil, Phys. Rev B 71, 125102 (2005). 
J. B. Kogut, Rev Mod. Phys. 51, 659 (1979). 

^' C. Dasgupta and B. Halperin, Phys. Rev Lett. 47, 1556 (1981). 
T. Senthil, A. Vishwanath, L. Balents, S. Sachdev, and M. P A. 
Fisher, Science 303, 1490 (2004). 

T. Senthil, L. Balents, S. Sachdev, A. Vishwanath, and M. P A. 
Fisher, Phys. Rev. B 70, 144407 (2004). 

S. Ki^agset, E. Sm0rgrav, J. Hove, F. S. Nogueira, and A. Sudb0, 
Phys. Rev Lett. 97, 247201 (2006). 

F. J. Jiang, M. Nyfeler, S. Chandrasekharan, and U. J. Wiese, J. 



16 



Stat. Mech. P02009 (2008). 

A. W. Sandvik, Phys. Rev. Lett. 98, 227202 (2007). 
" R. G. Melko and R. K. Kaul, Phys. Rev. Lett. 100, 017203 (2008). 
'^^ O. I. Motrunich and A. Vishwanath, arXiv:0805.1494 

A. B. Kuklov, M. Matsumoto, N. V. Prokof'ev, B. V. Svistunov, 

and M. Troyer, Phys. Rev. Lett. 101, 050405 (2008). 

D. Charrier, F. Alet, and P Pujol, Phys. Rev. Lett. 101, 167205 

(2008). 

^' S. Powell and J. T. Chalker, Phys. Rev. Lett. 101, 155702 (2008); 
Phys. Rev. B 78, 024422 (2008). 

A. F. Albuquerque et al., J. of Magn. and Magn. Materials 310, 
1187 (2007). 

M. Troyer, B. Ammon and E. Heeb, Lect. Notes Comput. Sci., 
1505, 191 (1998). 



